package scu.maqiang.numeric;

import scu.maqiang.mesh.ScalarFunc;

import java.util.Arrays;

public class MonteCarlo {
    public static void computePI() {
        int N = 80000000;
        int sum = 0;
        double x, y;
        for(int i = 0; i < N; i++){
            x = Math.random();
            y = Math.random();
            if( (x*x + y*y) < 1)
                sum++;
        }
        System.out.println("PI= " + (double)4*sum/N);

    }

    public static double MCIntegration1D(ScalarFunc f, double[] y, double[] domain, int N) {
        double length = Math.abs(domain[1] - domain[0]);
        double aveF = 0.0;
        double[] t = {0.0};
        for(int i = 0; i < N; i++) {
            t[0] = domain[0] + length * Math.random();
            aveF += f.action(y, 0, t);
        }
        return length * aveF / N;
    }

    public static Pair<double[], double[]> ODESolver(double[] inter, ScalarFunc func, double y0, int n) {
        double[] t = new double[n+1];
        double[] y = new double[n+1];
        double h = (inter[1] - inter[0]) / n;
        t[0] = inter[0];
        y[0] = y0;
        double[] time = {0};
        double[] yc = {0};
        double[] currentY = {0};
        for(int i = 0; i < n; i++) {
            time[0] = t[i];
            yc[0] = y[i];
            t[i+1] = t[i] + h;
            currentY[0] = y[i];
            y[i+1] = y[i] + MCIntegration1D(func, currentY, new double[]{t[i+1], t[i]}, 100000);
        }
        return new Pair<>(t, y);
    }

    public static void main(String[] args) {
        //computePI();
        ScalarFunc func = (xy, label, t) -> t[0] * t[0];
        double result = MCIntegration1D(func, null, new double[]{0, 2}, 100000);
        System.out.println("result = " + result);

        ScalarFunc f = (y, label, t) -> -1000 * y[0] + 3000 - 2000 * Math.exp(-t[0]);
        Pair<double[], double[]> result2 = ODESolver(new double[] {0, 1}, f, 0.0, 1000);
        System.out.println(Arrays.toString(result2.getFirst()));
        System.out.println(Arrays.toString(result2.getSecond()));
    }
}
